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Abstract. Linear response calculations based on the time-dependent density-functional theory are 
presented. Especially, we report results of the finite amplitude method which we have recently 
proposed as an alternative and feasible approach to the (quasiparticle-)random-phase approximation. 
Calculated properties of the giant resonances and low-energy El modes are discussed. We found a 
universal linear correlation between the low-energy E 1 strength and the neutron skin thickness. 
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INTRODUCTION 

Nuclei are created in stars via nuclear reactions where the temperatures are extremely 
high. Understanding of astrophysical phenomena generally requires many subfields of 
physics. Among them, nuclear physics plays key roles in the generation of elements, the 
evolution of the stars, the energy production in the universe, etc. The lightest elements up 
to helium were produced in the Big Bang. The other heavier elements are generated by 
many kinds of nuclear reactions in stars. Especially, in stellar explosions, there are thou- 
sands of reactions supposed to take place producing a variety of radioactive isotopes. 
However, to date, their reaction rates have not been hardly determined experimentally. 

To overcome these experimental difficulties, reliable theoretical information is nec- 
essary. Traditional ab-initio methods start from a nucleon-nucleon potential which de- 
scribes nucleon-nucleon scattering data. However, since the nuclear systems are strongly 
correlated because of a repulsive core in the potential, their description requires highly 
sophisticated many -body methods, such as the quantum monte carlo method [1]. To 
describe the nucleus in a quantitative way, they must employ an additional three-body 
force. These ab-initio methods are so involved that, even at present, their investigations 
have been limited to very light nuclei and to the homogeneous nuclear matter [2]. 

Under these circumstances, it is highly demanded to establish a universal theoretical 
approach which is able to describe properties of all species of nucleus. The nuclear 
density functional theory [3] is the most promising candidate among many nuclear 
models. It was referred to as the self-consistent mean-field model, but its concept 
is analogous to the density functional theory for many-electron systems [4, 5]. The 
functionals have typically about ten parameters which are adjusted by extensive fits to 



nuclear structure data. The most prominent feature of the model is that a single energy 
functional enables us to quantitatively describe almost all nuclei in the nuclear chart and 
infinite nuclear/neutron matter as well. 

The nuclear density functional models have been extensively used since 1970s [6]. In 
the beginning, there were several restrictions related to symmetries of the wave func- 
tions, which limited applications of the model. In 1990s, significant computational ad- 
vances together with vast amount of new spectroscopic data obtained with large gamma- 
ray arrays changed the situation. Using the cranking prescription, the density functional 
models were becoming a standard tool to study rotational bands in heavy nuclei [7]. In 
particular, the models were very successful in studies of superdeformed rotational bands 
at high spin in A = 150 [7, 8] and A = 190 regions [9]. It should be emphasized that the 
model parameters were never adjusted to these bands at all. In 2000s, systematic calcula- 
tions of nuclear ground- state properties were performed, to predict properties of nuclei 
far from the stability line [10]. At the same time, experimental developments greatly 
increased our knowledge on radioactive isotopes, then, had an impact on the energy 
density functionals. The nuclear density functional methods has now reached a point 
where one need to introduce "correlations beyond the Kohn-Sham scheme" to further 
improve the quality of the description. 

The density functional theory is designed for the description of ground- state prop- 
erties. Its extension to a time-dependent density functional theory is formally straight- 
forward, which is also analogous to the time-dependent mean-field theory. The density 
matrix p(t) obeys the following equation, 

ihj t p{t) = [h{t\p{t% (l) 

where h{t) is the Kohn-Sham (mean-field) Hamiltonian which is a functional of pit). 
In the case of open-shell superfluid nuclei with pairing, the density p (t) and the Hamil- 
tonian h(t) should be generalized to those with double size of matrices, R(t) and H(t), 
respectively [11]. Namely, R(t) contains not only the normal density p{t) but also the 
pairing tensor K(t). Similarly, H(t) has the pair potential A(t) in addition to h(t). It 
is not yet fully understood how we can compute excitation properties using the time- 
dependent density-functional theory, except for a few limiting cases. One such case is 
the low-energy regime of the nuclear dynamics, such as surface vibrations and shape 
fluctuations, which can be reached by the adiabatic limit of the time-dependent density- 
functional theory. It is closely related to the microscopic derivation of the Bohr model 
[12]. Recently, we have performed the numerical calculations of the large- amplitude 
quadrupole dynamics in various isotopes [13, 14, 15, 16]. 

Another important limiting case is the small-amplitude limit which provides us with 
a powerful tool to study linear response in nuclei. It is known as the random-phase ap- 
proximation (RPA) or the quasiparticle-random-phase approximation (QRPA) in nuclear 
physics, which accesses the regime of giant resonances. However, since the calculation 
of the residual interaction is tedious for the realistic energy functionals, it has been com- 
mon to ignore some terms in practice and to sacrifice the full self-consistency. To over- 
come these difficulties and facilitate an implementation of the full self-consistency, we 
employ two different methodologies; one is based on the real-time method (RTM) [17] 
and the other is the finite amplitude method (FAM) [18]. The RTM has an advantage 



that it does not require the calculation of the complex residual interactions, although 
it has a limitation in the achieved energy resolution inversely proportional to the time 
duration T. Our recent applications of the RTM are based on the canonical-basis frame- 
work [19] which is able to take into account dynamical pairing effects in nuclei. This 
is presented in another contribution to this volume [20]. In contrast, the FAM, which is 
complementary to the RTM, is a method of calculating the matrix elements of the resid- 
ual field, 8h = 8h/8p ■ 8p, using the finite difference. This does not require excessive 
programming but can be done by employing the program of the static density-functional 
calculation. We have performed systematic and fully self-consistent RPA calculations of 
photoabsorption cross sections for wide mass region (A < 100), for both spherical and 
deformed nuclei. In this report, we show results of these symmetry-unrestricted FAM 
calculations. 

FINITE AMPLITUDE METHOD 

In this section, we recapitulate the FAM. 



FAM for RPA 

First, we discuss the case without the pairing correlations. In this case, the FAM leads 
to residual fields appearing in the RPA. For more details, readers are referred to Ref . [18]. 

The linear-response RPA equation to a weak external field with a fixed frequency, 
Vext(w), can be expressed in terms of the forward and backward amplitudes, \Xj(o))) 
md (Yi(co)\. 

CO \Xi(G))) = (h - £i) \Xi((0)) +P{V cxt ((D) + 8h((0)} \<j>i), (2) 
-co (Yi(co)\ = (Yi(co)\ (h - £i) + {VextM + 8h((0)}P, (3) 

where the subscript i indicates the occupied orbitals (i = 1 , 2, • • • , A) and the operator P 
denotes the projector onto the particle space, P = 1 — Usually, the induced 

residual field 8h((o) is expanded to the first order with respect to \Xi((o)) and \Yi((o)). 
This leads to the well-known matrix form of the linear-response equation and calculation 
of these matrix elements is most time-consuming in practice. Instead, we utilize the fact 
that the linearization is numerically achieved for 8h((o) = h[po + 8p(co)] — ho within the 
linear approximation. In order to perform this numerical differentiation in the program, 
we use a small trick in the calculation of the single-particle Hamiltonian h[p]. 

First, we should notice that the Sh(oj) depends only on the forward "ket" amplitudes 
\Xi(co)) and backward "bra" ones (Yi(co)\. In other words, it is independent of bras 
(Xi(co)\ and kets \Yi(co)). This is related to the fact that the transition density 8p(co) 
depends only on \Xi(oj)} and (Yi(oj)\. 

Sp((D) = '£{\X l ((D))(<l> l \ + \<l> l )(Y l (ao)\}. (4) 



Then, we can calculate the residual fields in a following manner [18]: 

8h(G>) = ^(h\pn]-ho), (5) 

where 77 is a small real parameter to realize the linear approximation, p^ are defined by 
Pn^im + llM^m + riiYiia)])}. (6) 

i 

Once \Xi(co)) and (Yi(co) | are given, the calculation of h[pr,\ is an easy task. This does not 
require complicated programming, but only needs a small modification in the calculation 
of h[p] . Of course, eventually, we need to solve Eqs. (2) and (3) to determine the forward 
and backward amplitudes. We use an iterative algorithm to solve this problem. Namely, 

we start from initial amplitudes \X^) and then update them in every iteration, 

(pQ (w) ),(y/ n) |) ->• (\xj- n+1) ), (^ (w+1) |), until the convergence. In each step, we calculate 
8h{(o) using the FAM as Eq. (5). 



FAM for QRPA 

The FAM in the previous section can be extended to superfluid nuclei, namely, to the 
QRPA with the Bogoliubov extension of the mean fields. For more details, readers are 
referred to Ref. [21]. 

A self-consistent solution of static Kohn-Sham-Bogoliubov problems determines the 
ground-state densities (po, Kb) and the ground-state Hamiltonians (/z ,A ). They are 
given in terms of the quasiparticle wave functions, (l/^V^). Then, following the same 
argument in the previous section, we can derive equations for the residual fields, Sh(oo) 
and 8A(a) as follows [21]: 

8h{(o) = ^-{h[p rl ,K r} }-ho), 

" (V) 
5A ( W ) = -( a [Pt ?5 k 7? ]-A ), 

where the density and pairing tensor (p^, K>j) are defined by 

p n = (V* + i 1 UX)(V + riU*Y) T , (8) 
Kj] = (V* + riUX)(U + riV*Y) T . (9) 

Here, the forward and backward amplitudes (^ v ,y MV ) have subscripts jiv to specify 
two-quasiparticles. On the other hand, the subscripts of (U^, V^v) indicate a basis of the 
single-particle space (k) and the quasiparticle (/i). Again, utilizing an iterative algorithm 
for solution of the QRPA equation, we can solve the QRPA linear-response equation 
without explicitly calculating the residual interactions. 



Numerical results 



Development ofFAM computer programs 

We have developed the computer codes of the FAM in several representations. The 
FAM-RPA is available in the three-dimensional (3D) grid representation [22, 23]. This 
provides a completely symmetry unrestricted RPA calculation. All the single-particle 
wave functions and RPA amplitudes are represented by these grid points: 

Mr k , o)Mco;A, o)Ji{(o-r k , o)}[l\::: A ^= u ^ down . (10) 

Since the results are not sensitive to mesh spacing in a region outside of the interacting 
region, the adaptive grid representation is used to reduce the number of grid points Ng^ 
[17]. In the followings, we show results obtained with this code. 

The code for FAM-QRPA was also developed, but still with some symmetry restric- 
tions at present. In Ref. [21], we developed a code in the radial coordinate representation 
based on the spherical static program hfbrad [24], which assumes the spherical sym- 
metry in the ground state. Lately, another FAM-QRPA code has been developed in the 
harmonic oscillator representation [25], which is based on the program hfbtho with 
the axial symmetry [26]. 

In the present applications, we use complex energy, co = E + zT/2, with T = 1 
MeV, which introduces an artificial damping width. This smoothing is necessary in 
two reasons: To obtain smooth strength functions, and to speed up the convergence for 
the iterative procedure. However, the FAM formulae, Eqs. (5) and (7) themselves, are 
independent from the smoothing parameter T. In fact, we can use the FAM for explicit 
construction of the (Q)RPA matrix and calculate the discrete normal modes. This will 
be reported elsewhere [27]. 



Giant dipole resonances 

We have carried out the systematic calculation of electric dipole response in the FAM- 
RPA with the 3D grid representation. We show evolution of peak energies of the giant 
dipole resonances (GDR) as functions of mass number, neutron number, and proton 
number in Fig. 1. The GDR peak position is estimated by the average energy 

£peak= mi i C ° maX ^ m,( Wmax )= 1^ dCOO) k S{(0-D El ) (11) 

where S(co;Dei) = L„ IH-D£i|0)| 2 <5((0 — E n ). The maximum energy is (Omax ~ 40 
MeV. The El operator is defined with the recoil charges for protons, Ne/A and for 
neutrons, —Ze/A. In deformed nuclei, since the peak energy depends on direction of 
the El operator (x, y, and z), their averaged value is adopted in Fig. 1. The upper panel 
shows the GDR peak energies from oxygen to nickel, as a function of mass number. 
In the medium-mass region, the peak energies approximately follow the empirical low, 
21A -1 / 3 + 31A -1 / 6 MeV, denoted by the solid curve. However, in each isotopic chain, 
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FIGURE 1. Calculated GDR peak energies as functions of mass number (top), neutron number (bottom 
left), and proton number (bottom right). The isotopic chains are connected by lines for the top and bottom 
left panels, while the isotonic chains are shown in the bottom right. 
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FIGURE 2. Calculated PDR strength fraction as a function of neutron number. 



the peak energies in stable nuclei are the highest, while they are decreasing as leaving 
from the stability line. In addition, we can see some kind of shell effects. This can 
be more clearly seen in the lower panels of Fig. 1, in which the peak energies are 
plotted as functions of neutron and proton numbers. There are cusps at N =14 and 28 
corresponding to the subshell closure of ld 5 / 2 and I/7/2 orbitals. This may be attributed 
to the emergence of low-energy pygmy peaks beyond these neutron numbers (see the 
discussion in the next section). The proton shell effects seem not to be as significant as 
those of neutrons (see the bottom right panel of Fig. 1). 




Low-energy El strength 

Now, let us move to the low-energy part of the El strength distribution. In some 
nuclei, there appear small peaks in the El strength distribution, well separated from the 
main GDR peak, which are often called "pygmy dipole resonance" (PDR). In contrast 
to the main peak of GDR, the PDR strength distribution is sensitive to nuclear properties 
at nuclear surface and at low density. Thus, its property may provide us with useful 
constraints on the energy density functional, to identify the equation of state (EOS) of 
the nuclear and neutron matters. For instance, the neutron skin thickness is known to be 
well correlated with the slope of the neutron-matter EOS [28]. Thus, if the neutron skin 
thickness has a strong correlation with the low-energy El strength, we may pin down 
the EOS property by observing the PDR in experiments. 

First, we define the PDR strength fraction as 

/pdr = mi((O c )/mi(oo), (12) 

where m\((o) is given in Eq. (11) and we adopt co c = 10 MeV. In Fig. 2, we show the 
neutron-number dependence of /pdr- It indicates a strong shell effect. Namely, there 
are clear kinks at N = 14, 28, and 50. Let us concentrate our discussion on the kinks at 
N = 28. The PDR fractions suddenly increase at N = 28 — > 30 and continue to increase 
till N = 34 where the neutron 2p shell are filled. Beyond N = 34, the PDR fractions 
are roughly constant for 34 < N < 50, in which the neutrons are filling high-£ orbits of 
/5 /2 and gg n- Beyond N = 50, the neutrons start filling 2<i 5 / 2 orbits, then the /pdr again 
shows a sudden increase. These behaviors strongly suggest that the spatially extended 
nature of the low-£ neutron orbits near the Fermi level plays a primary role for the 
emergence and growth of the PDR. We have also observed that the deformation tends to 
increase the PDR strength, especially in the region N > 56. More detailed analysis can 
be found in our recent paper [23]. 

Finally, let us examine the correlation between /pdr and the neutron skin thickness. 
The skin thickness is defined by the difference in radius between neutrons and protons. 
Plotting the PDR fraction as a function of the skin thickness, we observe a linear 
correlation between them, but only in specific regions of the neutron number. This is 



illustrated in Fig. 3 for isotopes with Z = 16 22 and with Z = 26 32, which show 

the kinks at N = 28 and 50. The PDR fraction in each isotopic chain shows a linear 
correlation with the skin thickness in the regions of the neutron number N = 28 — 34 and 
N > 50. The positions of the kinks are located at different values of the skin thickness 
for different isotopes. However, the slope is universal for all the isotopes; 0.18 ~ 0.20 
fm _1 . Despite the fact that the deformation and shell ordering are different and vary 
from nucleus to nucleus, the universal linear correlation remains valid for 50 < N < 76. 
It should be noted that the linear correlation can be observed only for each isotopic 
chain. Deleting the lines connecting isotopic chains in Fig. 3, we only see scattered 
points showing a weak correlation. Again, for detailed analysis on this issue, readers are 
referred to Ref. [23]. 
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